Εργαστηριακή Άσκηση 2#

Σκοπός της δεύτερης σειράς ασκήσεων είναι η εξοικείωση με τις συναρτήσεις σχεδιασμού φίλτρων πεπερασμένης κρουστικής απόκρισης (FIR) και την υλοποίησή τους στο MATLAB. Προτού ξεκινήσετε την άσκηση θα πρέπει να μελετήσετε με προσοχή το Κεφάλαιο 1 και, ειδικότερα, την παράγραφο 1.3 του τεύχους του μαθήματος. Το MATLAB (www.mathworks.com) είναι ένα διαδραστικό εμπορικό πρόγραμμα (Windows, Linux, Unix) με το οποίο μπορείτε να κάνετε εύκολα αριθμητικές πράξεις με πίνακες. Μπορεί να το έχετε εγκατεστημένο τοπικά, στον προσωπικό σας υπολογιστή ή να εργάζεστε σε κάποιο Εργαστήριο Προσωπικών Υπολογιστών (ΕΠΥ) της Σχολής σας που διαθέτει το συγκεκριμένο λογισμικό.

Live Code

Press the following button to make python code interactive. It will connect you to a kernel once it says “ready” (might take a bit, especially the first time it runs).

Μέρος 1: Εισαγωγή#

Στο MATLAB οι συναρτήσεις \(fft\) και \(ifft\) υποθέτουν ζεύγος μετασχηματισμού Fourier \(x(t)\) και \(X(f)\) υπολογισμένων σε μη αρνητικά διαστήματα \( t=[0:N-1]ts \) και \( f=[0:N-1]fo. \) Όπως έχετε ήδη δει στην Εργαστηριακή άσκηση 1, το άνω μισό μέρος του διαστήματος συχνοτήτων αντιστοιχεί στις αρνητικές συχνότητες του σήματος, όταν υπολογίζουμε το \(X(f)\) με τη βοήθεια της συνάρτησης \(fft\) (για σήματα πραγματικών τιμών, αυτά τα δυο μισά είναι κατοπτρικά ως προς το μέσο του διαστήματος). Ακριβώς το ίδιο ισχύει και για το άνω μισό μέρος του χρονικού διαστήματος, όταν το σήμα \(x(t)\) προκύπτει από τον αντίστροφο μετασχηματισμό Fourier μέσω της \(ifft\).

Το MATLAB διαθέτει τη συνάρτηση \(fftshift\) για να ολισθήσει κυκλικά τις τιμές του σήματος ή του μετασχηματισμού Fourier, ώστε να αντιστοιχούν σε κεντραρισμένα στο μηδέν αμφίπλευρα διαστήματα, δηλαδή, στις χρονικές στιγμές \(tb=[-ceil((N-1)/2): floor((N-1)/2)]ts\) ή στις συχνότητες \(fb=[–ceil((N-1)/2): floor((N-1)/2)]fo\). Με τον τρόπο αυτό μπορούμε να παράγουμε τα \(xb(t)\) και \(Xb(f)\) που αντιστοιχούν στην αμφίπλευρη αναπαράσταση του σήματος και του μετασχηματισμού Fourier.

Για να κατανοήσετε τα ανωτέρω θεωρείστε το διάνυσμα \([1 2 3 4]\) ως το αποτέλεσμα του \(FFT\) μήκους 4. Τότε, το πρώτο στοιχείο (1) είναι ο όρος dc, το τρίτο στοιχείο (3) είναι το σημείο στο μισό της συχνότητας δειγματοληψίας \(fs/2\), που μπορεί να εκληφθεί ότι αντιστοιχεί είτε στην \(–fs/2\) είτε στην \(+fs/2\). Τα στοιχεία 2 και 4 αντιστοιχούν στις συχνότητες \(+fs/4\) και \(–fs/4\). Εφαρμόζοντας την \(fftshift\), το στοιχείο 3 εμφανίζεται πρώτο, που σημαίνει ότι στο MATLAB αντιστοιχεί στην αρνητική συχνότητα \(–fs/2\), το επόμενο στοιχείο 4 αντιστοιχεί στη συχνότητα \(–fs/4\) ακολουθούμενο από το dc και τη συχνότητα \(+fs/4\). Για ένα μετασχηματισμό περιττού μήκους, δεν υφίσταται σημείο για το \(±fs/2\). Έτσι για το διάνυσμα \([1 2 3]\), η εφαρμογή της \(fftshift\) θα δώσει τα στοιχεία που αντιστοιχούν στις συχνότητες \(–fs/3, 0, +fs/3\).

Εκτός του ότι παράγουν εξόδους με τις αρνητικές συχνότητες ή χρόνους στο άνω μισό του διανύσματος, αμφότερες οι συναρτήσεις \(fft\) και \(ifft\) αναμένουν ως είσοδο διάνυσμα με την ίδια μορφή, αφού προφανώς ισχύουν οι ταυτότητες

\(h == \text{ifft}(\text{fft}(h))\) και \(H == \text{fft}(\text{ifft}(H))\)

Η πρώτη υποδεικνύει ότι η είσοδος της \(ifft\) πρέπει να είναι αντεστραμμένη, όπως την παράγει η \(fft\), και η δεύτερη ότι η είσοδος της \(fft\) πρέπει να είναι αντεστραμμένη, όπως την παράγει η \(ifft\).

Στο επόμενο σχήμα βλέπετε παραστατικά ένα ημιτονικό σήμα που έχει πολλαπλασιασθεί με παράθυρο Blackman τόσο στην αμφίπλευρη, όσο και την μονόπλευρη αναπαράστασή του.

import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets
from ipywidgets import Layout
from IPython.display import display, clear_output

np.random.seed(0)

# Generate random data
t = np.linspace(-40, 40, 80)  # 80 samples
random_data = np.random.randn(len(t))

# Initial start and end points for the window
start, end = 30, 50

# Output widget for the plots
output = widgets.Output()

def apply_window_and_plot(start, end):
    with output:
        clear_output(wait=True)
        # Create a sinusoidal window
        sin_window = np.zeros_like(t)
        sin_window[start:end] = np.sin(np.linspace(0, np.pi, end - start))

        # Apply sinusoidal window to random data
        windowed_data = sin_window * random_data

        # Append the negative times to the end of the time array
        t_append_negative = np.concatenate((t[t >= 0], t[t < 0] + t[-1] - t[0] + (t[1] - t[0])))
        windowed_data_append_negative = np.concatenate((windowed_data[t >= 0], windowed_data[t < 0]))
        plt.close('all')  # Close all existing figures
        # Create subplots
        fig, axs = plt.subplots(2, 1, figsize=(10, 10))

        # Plot original windowed data
        axs[0].stem(t, windowed_data, linefmt='C0-', markerfmt='C0o', basefmt=" ")
        axs[0].set_title("Original Sinusoidal Windowed Random Data", fontsize=16)
        axs[0].set_xlabel("Time (samples)", fontsize=12)
        axs[0].set_ylabel("Amplitude", fontsize=12)

        # Plot with negative times appended
        axs[1].stem(t_append_negative, windowed_data_append_negative, linefmt='C1-', markerfmt='C1o', basefmt=" ")
        axs[1].set_title("Sinusoidal Windowed Random Data with Negative Times Appended", fontsize=16)
        axs[1].set_xlabel("Time (samples)", fontsize=12)
        axs[1].set_ylabel("Amplitude", fontsize=12)

        # Adjust layout and show plot
        plt.tight_layout()
        
        plt.show()
        

# Define the range slider
window_range_slider = widgets.IntRangeSlider(
    value=[start, end],
    min=0,
    max=len(t),
    step=1,
    description='Window Range:',
    layout=Layout(width='90%'),  # Setting the width to 90%
    continuous_update=False
)

# Function to update the plot based on the range slider
def on_range_change(change):
    new_start, new_end = change['new']
    apply_window_and_plot(new_start, new_end)

# Observe the range slider for changes
window_range_slider.observe(on_range_change, names='value')

# Group the slider and output widget in a vertical box layout
vbox = widgets.VBox([window_range_slider, output])

# Display the VBox
display(vbox)

# Initialize with the default plot
apply_window_and_plot(start, end)
<<<<<<< HEAD
=======
>>>>>>> d546939e94d695c6e9c44c085bc5ec2e2e14440e

Οι αντίστοιχοι μετασχηματισμοί Fourier, σε αμφίπλευρη και μονόπλευρη αναπαράσταση, φαίνονται στο επόμενο σχήμα:

lab2_2.png
../_images/e140b64bc6b77ff76c35282e68bbf0babbae27f2b55e4369b432a0cabea23679.png ../_images/c8fa559816efd87255cb0408803a32e0d9122c4c5477e8cf97b7ae457648bdc6.png

Όταν τα \(x(t)\) και \(X(f)\) παράγονται από το MATLAB, δεν χρειάζεται κάποια ιδιαίτερη προσοχή, πλην της κυκλικής ολίσθησης σε περίπτωση που θέλουμε π.χ. να σχεδιάσουμε το αμφίπλευρο φάσμα ή σήμα. Όταν όμως ένα εκ των \(x(t)\) ή \(X(f)\) ορίζεται από τον χρήστη απαιτείται περισσότερη προσοχή, διότι, συνήθως χρησιμοποιούνται τα αμφίπλευρα σήματα ή φάσματα. Μπορείτε να μεταβείτε από τη μία αναπαράσταση στην άλλη ως εξής:

\(x = ifftshift(xb)\), \(X = fft(x)\), \(Xb = fftshift(X)\), εάν ξεκινάτε από αμφίπλευρο σήμα και θέλετε να καταλήξετε σε αμφίπλευρο φάσμα, και

\(X = ifftshift(Xb)\), \(x = ifft(X)\), \(xb = fftshift(x)\), εάν ξεκινάτε από αμφίπλευρο φάσμα και θέλετε να καταλήξετε σε αμφίπλευρο σήμα,

όπου η συνάρτηση \(ifftshift\) του MATLAB εκτελεί την αντίστροφη λειτουργία της \(fftshift\). Όταν το \(N\) είναι άρτιο, οι \(fftshift\) και \(ifftshift\) δίνουν το ίδιο αποτέλεσμα. Όταν όμως το \(N\) είναι περιττό αυτό δεν ισχύει και χρειάζεται προσοχή στη χρήση τους. Στην πράξη, η προσεκτική εφαρμογή των ανωτέρω έχει σημασία όταν υπολογίζεται η φάση του φάσματος. Το πλάτος του φάσματος δεν επηρεάζεται από την κυκλική ολίσθηση των στοιχείων που προκαλούν οι \(fftshift\) και \(ifftshift\) (δείτε ιδιότητες DFT).

Εξάσκηση#

Δοκιμάστε στο παράθυρο εντολών τα ακόλουθα προκειμένου να εμπεδώσετε τη χρήση των συναρτήσεων \(fftshift\) και \(ifftshift\).

import numpy as np

# Create the array X
X = np.arange(-2, 3)  # equivalent to MATLAB's -2:2

# Apply fftshift and ifftshift
fftshifted_X = np.fft.fftshift(X)
ifftshifted_X = np.fft.ifftshift(X)

# Double fftshift and a combination of fftshift and ifftshift
Y = np.fft.fftshift(np.fft.fftshift(X))
Z = np.fft.ifftshift(np.fft.fftshift(X))

# Check if arrays are equal
X_equals_Y = np.array_equal(X, Y)
X_equals_Z = np.array_equal(X, Z)

print(f"X equals Y: {X_equals_Y}")
print(f"X equals Z: {X_equals_Z}")
X = [-2:2]
fftshift(X)
ifftshift(X)
Y = fftshift(fftshift(X));
Z = ifftshift(fftshift(X));
isequal(X,Y)
isequal(X,Z)
X equals Y: False
X equals Z: True

Ερώτηση 1: Ποιο εκ των διανυσμάτων Υ και Ζ ισούται με το X; Γράψτε την απάντησή σας σε ένα αρχείο κειμένου lab2_nnnnn.txt, όπου nnnnn τα πέντε τελευταία νούμερα του αριθμού μητρώου σας, χρησιμοποιώντας το Notepad από το μενού των Windows (Start → Programs → Accessories → Notepad) και αποθηκεύστε το στον φάκελο My Documents. Θα υποβάλετε το αρχείο αυτό ηλεκτρονικά στο τέλος, αφού απαντήσετε και τις επόμενες ερωτήσεις, οπότε μπορείτε να τα αφήσετε ανοικτό.

Ερώτηση 2: Επαναλάβατε με \(X=[-1:2]\). Τι παρατηρείτε; Γράψτε την απάντησή σας στο αρχείο κειμένου lab2_nnnnn.txt. Δοκιμάστε στο παράθυρο εντολών τα ακόλουθα δύο παραδείγματα για να εμπεδώσετε τη χρήση των συναρτήσεων \(fftshift\) και \(ifftshift\) σε συνδυασμό με τις \(fft\) και \(ifft\).

Code
import numpy as np
from scipy.fft import fft, ifft, fftshift, ifftshift
import matplotlib.pyplot as plt
from ipywidgets import interact, Layout
import ipywidgets as widgets

# First part: Original signal and its FFT
xb = np.array([1, 2, 3, 4, 5, 4, 3, 2, 1])
x = ifftshift(xb)
X = fft(x)
Xb = fftshift(X)  # Spectrum with DC component in the center

# Second part: Low-pass filter effect
Xb_low_pass = np.array([0, 0, 1, 1, 1, 1, 1, 0, 0])
X_low_pass = ifftshift(Xb_low_pass)
x_low_pass = ifft(X_low_pass)
xb_low_pass = fftshift(x_low_pass.real)

# Function to plot
def plot_signals():
    fig1, axs = plt.subplots(2, 2, figsize=(15, 10))

    time_range = np.arange(-4, 5)
    
    axs[0, 0].stem(time_range, xb, linefmt='blue', markerfmt='bo', basefmt='r-')
    axs[0, 0].set_title('Original Signal')
    axs[0, 0].set_xlabel('Time (s)')
    axs[0, 0].set_ylabel('Amplitude')

    axs[0, 1].stem(time_range, np.abs(Xb), linefmt='red', markerfmt='ro', basefmt='r-')
    axs[0, 1].set_title('Magnitude Spectrum')
    axs[0, 1].set_xlabel('Frequency (Hz)')
    axs[0, 1].set_ylabel('Magnitude')

    axs[1, 0].stem(time_range, Xb_low_pass, linefmt='green', markerfmt='go', basefmt='r-')
    axs[1, 0].set_title('Low-pass Spectrum')
    axs[1, 0].set_xlabel('Frequency (Hz)')
    axs[1, 0].set_ylabel('Magnitude')

    axs[1, 1].stem(time_range, xb_low_pass, linefmt='purple', markerfmt='mo', basefmt='r-')
    axs[1, 1].set_title('Reconstructed Signal')
    axs[1, 1].set_xlabel('Time (s)')
    axs[1, 1].set_ylabel('Amplitude')

    plt.tight_layout()
    plt.show()

plot_signals()
close all; clear all;clc;
xb=[1 2 3 4 5 4 3 2 1] % πραγματικό σήμα με άρτια συμμετρία
figure; subplot (2,1,1); plot([-4:4],xb); ylabel('xb');
x=ifftshift(xb) % το σήμα με τις αρνητικές συνιστώσες στο άνω μέρος
X=fft(x) % FFT
Xb=fftshift(X) % το φάσμα με τη dc συνιστώσα στο κέντρο, πραγματικές
% τιμές με άρτια συμμετρία όπως αναμένεται
subplot (2,1,2); plot([-4:4],Xb);ylabel('Xb');
close all; clear all;clc;
Xb=[0 0 1 1 1 1 1 0 0] % φάσμα βαθυπερατού σήματος με άρτια συμμετρία
figure; subplot (2,1,1); plot([-4:4],Xb); ylabel('Xb');
X=ifftshift(Xb) % το φάσμα με τις αρνητικές συνιστώσες στο άνω μέρος
x=ifft(X) % IFFT
xb=fftshift(x) % πραγματικό σήμα με άρτια συμμετρία όπως αναμένεται
subplot (2,1,2); plot([-4:4],xb); ylabel('xb');
../_images/0bfb4e092a51e5dc4be78516e9d5366aed3c2c516f247e99c842c9605ab0da74.png

Ερώτηση 3: Τροποποιείστε το προηγούμενο παράδειγμα ώστε να ξεκινήσετε απευθείας με τον ορισμό του φάσματος του βαθυπερατού σήματος \(X\) όπως το αναμένει η \(ifft\). Γράψτε την απάντησή σας στο αρχείο κειμένου lab2_nnnnn.txt.

Μέρος 2: Σχεδιασμός και υλοποίηση φίλτρων#

Θα ασχοληθούμε με το Παράδειγμα 1.2 της παραγράφου 1.5 του τεύχους Μαθήματος. Το παράδειγμα αυτό παρουσιάζει δύο εναλλακτικές μεθόδους σχεδιασμού FIR φίλτρων:

α) τη μέθοδο των παραθύρων και

β) τη μέθοδο των ισοϋψών κυματώσεων

τις οποίες εφαρμόζει στην περίπτωση βαθυπερατών φίλτρων.

Ο κώδικας σε Matlab περιέχεται στο βιβλίο του μαθήματος στο παράδειγμα 1.2 (Κώδικας 1.3). Στο παράδειγμα, τα φίλτρα δοκιμάζονται σε ένα πραγματικό σήμα, s, το οποίο είναι αποθηκευμένο στο αρχείο sima.mat (binary αρχείο MATLAB). Πρόκειται για ένα σήμα sonar με φάσμα που εκτείνεται μέχρι περίπου τα 4 KHz και συχνότητα δειγματοληψίας Fs=8192 (είναι και αυτή αποθηκευμένη στο αρχείο sima.mat, μαζί με το σήμα). Παρατηρήστε ότι, στον κώδικα του παραδείγματος, προηγείται η φόρτωση του sima.mat (γραμμή 6), γεγονός που επιτρέπει τη χρήση των μεταβλητών s, Fs στο υπόλοιπο τμήμα του (γραμμές 7-44).

Η μέθοδος των παραθύρων#

Η μέθοδος των παραθύρων εφαρμόζεται στις γραμμές 8-38 του κώδικα. Πρώτο βήμα αποτελεί ο ορισμός της απόκρισης συχνότητας H ενός ιδανικού βαθυπερατού φίλτρου με συχνότητα αποκοπής \(Fs/8\) (γραμμή 9), μέσω ενός διανύσματος μήκους \(Ν=Fs\), γεγονός που οδηγεί σε ανάλυση συχνότητας \(fο = Fs /N = 1 Hz\). Το διάνυσμα H αποτελείται από μια αλληλουχία μονάδων και μηδενικών που δημιουργούνται από την κλήση των συναρτήσεων \(ones\) και \(zeros\), αντίστοιχα. Για περισσότερες πληροφορίες σχετικά με αυτές τις συναρτήσεις μπορείτε να συμβουλευτείτε την τεκμηρίωση του MATLAB, πληκτρολογώντας \(doc ‘function-name’\) στο παράθυρο εντολών, όπου \(‘function-name’\) το όνομα της συνάρτησης. Τα \(Fs/8\) πρώτα στοιχεία του \(H\), που αντιστοιχούν στις συχνότητες \([0,Fs/8)\), είναι μονάδες, ακολουθούν \(3Fs/4\) μηδενικά και άλλες \(Fs/8\)\( μονάδες που αντιστοιχούν στη ζώνη συχνοτήτων \)[–Fs/8, 0)$.

Επόμενο βήμα είναι ο υπολογισμός του αντίστροφου διακριτού μετασχηματισμού \(Fourier\) \((IDFT)\), που υπολογίζεται από τη συνάρτηση \(ifft\) του MATLAB (γραμμή 12). Ακολουθεί μια αναδιάταξη του αποτελέσματος του αντίστροφου \(DFT\) (γραμμές 13-14), που ισοδυναμεί με ολίσθηση της κρουστικής απόκρισης του φίλτρου κατά το μισό της μήκος. Στη συνέχεια, η κρουστική απόκριση \(h\) περικόπτεται σε μήκη 32+1, 64+1 και 128+1 δειγμάτων (γραμμές 15-17). Με τη βοήθεια του \(wvtool\) (Window Visualization Tool), απεικονίζονται στο ίδιο διάγραμμα η κρουστική απόκριση και η απόκριση συχνότητας του βαθυπερατού φίλτρου και για τα 3 παραπάνω μήκη. Το \(wvtool\) είναι ένα παρόμοιο εργαλείο με το \(wintool\) και χρησιμεύει για την οπτικοποίηση παραθύρων στο πεδίο του χρόνου και της συχνότητας. Παρατηρήστε ότι όσο μεγαλύτερο το μήκος του φίλτρου τόσο μικρότεροι είναι οι πλευρικοί λοβοί στην απόκριση συχνότητας. Για να μειωθούν ακόμα περισσότερο αυτοί οι πλευρικοί λοβοί και οι επιπτώσεις του ορθογωνικού παραθύρου, κατασκευάζονται παράθυρα \(Hamming\) και \(Kaiser\) (μέσω των συναρτήσεων \(hamming\) και \(kaiser\) αντίστοιχα, γραμμές 24-25) τα οποία εφαρμόζονται στο φίλτρο μήκους 64+1 σημείων \(h64\) (γραμμές 27-30). Παρατηρήστε μέσω του \(wvtool\) (γραμμή 31) τη σαφώς χαμηλότερη στάθμη των πλευρικών λοβών στην απόκριση συχνότητας των παραθύρων \(Hamming\) και \(Kaiser\) σε σχέση με το ορθογωνικό. Τέλος, το σήμα s φιλτράρεται με καθένα από τα τρία φίλτρα (ορθογωνικό, \(Hamming\), \(Kaiser\)), χρησιμοποιώντας τη συνάρτηση \(conv()\) που υπολογίζει τη συνέλιξη μεταξύ του σήματος και της κρουστικής απόκρισης του εκάστοτε φίλτρου (γραμμές 33-38). Παράλληλα, υπολογίζεται και σχεδιάζεται η πυκνότητα φάσματος ισχύος κατά \(Welch\) με τη βοήθεια της συνάρτησης \(pwelch\). Το αποτέλεσμα δείχνει εμφανώς την ραγδαία μείωση της φασματικής ισχύος για συχνότητες άνω της \(Fs/8\).

Μέθοδος ισοϋψών κυματώσεων Η μέθοδος των ισοϋψών κυματώσεων εφαρμόζεται στις γραμμές 40-44 του κώδικα για την κατασκευή ενός φίλτρου \(Parks-McClellan\). Η κατασκευή του φίλτρου πραγματοποιείται με την κλήση της συνάρτησης \(firpm\) (γραμμή 41). Για τον ορισμό των παραμέτρων εισόδου της \(firpm\) συμβουλευτείτε την τεκμηρίωση του MATLAB. Η εφαρμογή του φίλτρου στο σήμα γίνεται όπως και στην προηγούμενη μέθοδο, με χρήση της συνάρτησης conv (γραμμή 43).

Code
import numpy as np
from scipy.signal import remez, convolve, welch
import matplotlib.pyplot as plt
import scipy.io
import sounddevice as sd

#6
with open('../files/sima.txt') as f:
    s = [float(x) for x in f]
s=np.array(s)
Fs=8192


# 7
f, Pxx = scipy.signal.welch(s, Fs)
plt.figure(); plt.semilogy(f, Pxx); plt.grid(True); plt.xlabel('Frequency [Hz]'); plt.ylabel('PSD [V**2/Hz]')

# 8-9
N = Fs // 2  # Adjusted for Python indexing, assuming Fs is even
H = np.concatenate((np.ones(N//4), np.zeros(N//2), np.ones(N//4)))

# 10-14
h = np.fft.ifft(H, n=N).real
h = np.fft.fftshift(h)  # Equivalent to MATLAB's circshift for centering the impulse response

# 15-17
h32 = h[N//2-16:N//2+17]
h64 = h[N//2-32:N//2+33]
h128 = h[N//2-64:N//2+65]

# 18 - Stem plot for h64
plt.figure()
plt.stem(range(len(h64)), h64, basefmt=" ")
plt.title('Stem Plot of h64')
plt.xlabel('Samples')
plt.ylabel('Amplitude')
plt.grid()

# 19 - Frequency response of h64
w, h = scipy.signal.freqz(h64)
plt.figure()
plt.plot(w / np.pi * (Fs / 2), np.abs(h))  # Normalizing frequency to Hz
plt.title('Frequency Response of h64')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.grid()

# 20 - For visualizing the frequency responses of h32, h64, and h128 together
plt.figure(figsize=(14, 4))

for filt, label in zip([h32, h64, h128], ['h32', 'h64', 'h128']):
    w, h = scipy.signal.freqz(filt)
    plt.plot(w / np.pi * (Fs / 2), 20 * np.log10(abs(h)), label=label)

plt.title('Frequency Responses of h32, h64, h128')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude (dB)')
plt.legend()
plt.grid()
plt.grid(True)
plt.show()

# 22-26
wh = np.hamming(len(h64))
wk = np.kaiser(len(h64), 5)
plt.figure(); plt.plot(wh, 'b', label='Hamming'); plt.plot(wk, 'r', label='Kaiser'); plt.legend(); plt.grid()

# 27-30
h_hamming = h64 * wh
h_kaiser = h64 * wk

# 31 `wvtool` equivalent visualization
# Plot the impulse responses
plt.figure(figsize=(10, 6))
markers, stemlines, baseline = plt.stem(range(len(h64)), h64, linefmt='C0-', markerfmt='C0o', basefmt="C0", label='h64')
plt.setp(markers, markersize = 5)  # Adjusting marker size for better visibility
plt.setp(stemlines, linestyle = '-')  # Adjusting stem line style

markers, stemlines, baseline = plt.stem(range(len(h_hamming)), h_hamming, linefmt='C1-', markerfmt='C1o', basefmt="C1", label='h_hamming')
plt.setp(markers, markersize = 5)
plt.setp(stemlines, linestyle = '-')

markers, stemlines, baseline = plt.stem(range(len(h_kaiser)), h_kaiser, linefmt='C2-', markerfmt='C2o', basefmt="C2", label='h_kaiser')
plt.setp(markers, markersize = 5)
plt.setp(stemlines, linestyle = '-')

plt.title('Impulse Responses')
plt.xlabel('Samples')
plt.ylabel('Amplitude')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# Plot the frequency responses
plt.figure(figsize=(14, 4))
for filt, label, color in zip([h64, h_hamming, h_kaiser], ['h64', 'h_hamming', 'h_kaiser'], ['C0', 'C1', 'C2']):
    w, h = scipy.signal.freqz(filt)
    plt.plot(w / np.pi * (Fs / 2), 20 * np.log10(abs(h)), label=label, color=color)

plt.title('Frequency Responses')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude (dB)')
plt.legend()
plt.grid()
plt.show()

# 32-38 Filtering and plotting PSD
y_rect = np.convolve(s, h64, mode='same')
f, Pxx = scipy.signal.welch(y_rect, Fs)
plt.figure();plt.grid(True); plt.semilogy(f, Pxx); plt.xlabel('Frequency [Hz]'); plt.ylabel('PSD [V**2/Hz]'); plt.title('h64')

y_hamm = np.convolve(s, h_hamming, mode='same')
f, Pxx = scipy.signal.welch(y_hamm, Fs)
plt.figure();plt.grid(True); plt.semilogy(f, Pxx); plt.xlabel('Frequency [Hz]'); plt.ylabel('PSD [V**2/Hz]'); plt.title('h_hamming')

y_kais = np.convolve(s, h_kaiser, mode='same')
f, Pxx = scipy.signal.welch(y_kais, Fs)
plt.figure(); plt.grid(True); plt.semilogy(f, Pxx); plt.xlabel('Frequency [Hz]'); plt.ylabel('PSD [V**2/Hz]'); plt.title('h_kaiser')

# 40-44  Parks-McClellan low-pass filter
freq = np.array([0, 0.10, 0.15, 0.5]) * Fs  # Scale frequencies correctly
gain = [1, 1, 0, 0]

# Call firwin2 with corrected parameters and using the `fs` parameter
hpm = signal.firwin2(65, freq, gain, fs=Fs)
s_pm = np.convolve(s, hpm, 'same')

# Plotting or analyzing the filtered signals (y_rect, y_hamm, y_kais, s_pm) can be done similarly
# For example, to plot the PSD:
f, Pxx_spec = signal.welch(s_pm, Fs, 'flattop', 1024, scaling='spectrum')
plt.figure()
plt.semilogy(f, np.sqrt(Pxx_spec))
plt.xlabel('frequency [Hz]')
plt.ylabel('Linear spectrum [V RMS]')
plt.title('Power spectrum (welch)')
plt.show()

# 45-46 Playing sound (s and s_pm assumed as original and filtered signal)
sd.play(20 * s, Fs); sd.wait()
# sd.play(20 * s_lp, Fs); sd.wait()  # Uncomment and adjust variable name as needed for the filtered signal
1. clear all; close all;
2. % Το αρχείο "sima.mat" περιέχει το σήμα s και τη συχνότητα
3. % δειγματοληψίας Fs. Το φάσμα του σήματος εκτείνεται σχεδόν σε όλη την
4. % περιοχή συχνοτήτων μέχρι 4 KHz. Πάνω από 1 KHz, όμως, είναι θόρυβος
5. % και πρέπει να φιλτραριστεί.
6. load sima;
7. figure; pwelch(s,[],[],[],Fs);
8. % Ορίζεται η ιδανική βαθυπερατή συνάρτηση Η, με συχνότ. αποκοπ. Fs/8
9. H=[ones(1,Fs/8) zeros(1,Fs-Fs/4) ones(1,Fs/8)];
10. % Υπολογίζεται η κρουστική απόκριση με αντίστροφο μετασχ. Fourier
11. % Εναλλακτικά, μπορεί να χρησιμοποιηθεί η αναλυτική σχέση Sa(x)
12. h=ifft(H,'symmetric');
13. middle=length(h)/2;
14. h=[h(middle+1:end) h(1:middle)];
15. h32=h(middle+1-16:middle+17);
16. h64=h(middle+1-32:middle+33);
17. h128=h(middle+1-64:middle+65);
18. % figure; stem([0:length(h64)-1],h64); grid;
19. % figure; freqz(h64,1); % σχεδιάζουμε την απόκριση συχνότητας της h64
20. wvtool(h32,h64,h128); % αποκρίσεις συχνότητας των περικομμένων h
21. % Οι πλευρικοί λοβοί είναι υψηλοί!
22. % Πολλαπλασιάζουμε την περικομμένη κρουστική απόκριση με κατάλληλο
23. % παράθυρο. Χρησιμοποιούμε την h64 και παράθυρα hamming και kaiser
24. wh=hamming(length(h64));
25. wk=kaiser(length(h64),5);
26. figure; plot(0:64,wk,'r',0:64,wh,'b'); grid;
27. h_hamming=h64.*wh';
28. % figure; stem([0:length(h64)-1],h_hamming); grid;
29. % figure; freqz(h_hamming,1);
30. h_kaiser=h64.*wk';
31. wvtool(h64,h_hamming,h_kaiser);
32. % Φιλτράρουμε το σήμα μας με καθένα από τα τρία φίλτρα
33. y_rect=conv(s,h64);
34. figure; pwelch(y_rect,[],[],[],Fs);
35. y_hamm=conv(s,h_hamming);
36. figure; pwelch(y_hamm,[],[],[],Fs);
37. y_kais=conv(s,h_kaiser);
38. figure; pwelch(y_kais,[],[],[],Fs);
39. %
40. % Βαθυπερατό Parks-MacClellan
41. hpm=firpm(64, [0 0.10 0.15 0.5]*2, [1 1 0 0]);
42. % figure; freqz(hpm,1);
43. s_pm=conv(s,hpm);
44. figure; pwelch(s_pm,[],[],[],Fs);
45. sound(20*s); % ακούμε το αρχικό σήμα, s
46. sound(20*s_lp); % ακούμε το φιλτραρισμένο σήμα, s_lp
../_images/19c9114436945a38713023d244136939a41626ebfa3ec4ea6818c415b507402e.png ../_images/e1524a529ed04ac59de2baa268e26ea83927b092eb750ca469ae4c426b1c8f43.png ../_images/058b252f161df67cc8cc238a36d99eb3586acf6dd0c6a463c2c44dc84cc62dff.png ../_images/9bddee962a55b230f0ea49d59f88160cb4e64dab8b4091c06c9a7d575f61eece.png ../_images/4a72d6151a36a30c97a17200b0f47d473c34f5fba55c201ecec23b1ce1a3c709.png ../_images/44bcbdb88d45d63aa0a874b19c2ad715689a6041a78db9fe9e92dee44e8c561b.png ../_images/30c988880f69d37a68353a8d9645028e43c09337ea2c5995e357c2e6b97c822a.png ../_images/d7808dedd65e7989b17a62d043aa03539e98c23297e0170a7bd4b8fc83a5daa2.png ../_images/2fad1bc832295348cf429f103441705e4351e7524a9c934fe1c0f7644d3efb35.png ../_images/33bac07ce094128d53e69c9fac4148566bf0f0b5b9006f46972f822fc358bf6c.png ../_images/556eb9c777cda89af346e9ce6b50015d84c91a03db49952abfc0cb6ac2da42a1.png

Πειραματισθείτε#

  1. Τροποποιείστε τον κώδικα στην γραμμή 14 ώστε να επιτύχετε το ίδιο αποτέλεσμα με τη χρήση μιας εκ των συναρτήσεων ifftshift ή fftshift.

  1. Τροποποιείστε τον κώδικα ώστε να χρησιμοποιηθεί βαθυπερατό φίλτρο μήκους 140+1 αντί του 64+1. Σχεδιάστε την απόκριση συχνότητας του φίλτρου με σχεδιασμό παραθύρου, όπως στο παράδειγμα, για δύο περιπτώσεις: ορθογωνικού και Hamming.

Ορθογωνικό παράθυρο (απλή περικοπή της h)#

%matplotlib inline
import ipywidgets as widgets
from IPython.display import display, clear_output
import matplotlib.pyplot as plt
import numpy as np

# Assuming H is defined and calculated as before
H = np.hstack((np.ones(int(Fs/8)), np.zeros(int(Fs-Fs/4)), np.ones(int(Fs/8))))
h = np.real(np.fft.ifft(H))
middle = int(len(h)/2)
h = np.hstack((h[middle:], h[:middle]))
h32 = h[middle-16:middle+16]
h64 = h[middle-32:middle+32]
h128 = h[middle-64:middle+64]
h160 = h[middle-80:middle+80]
h256 = h[middle-128:middle+128]

h_variants = {
    'h32': h[middle-16:middle+16],
    'h64': h[middle-32:middle+32],
    'h128': h[middle-64:middle+64],
    'h160': h[middle-80:middle+80],
    'h256': h[middle-128:middle+128],
}

output2 = widgets.Output()

# Function to plot
def plot_stem(h_key):
    with output2:
        clear_output(wait=True)  # Clear the previous plots
        h_data = h_variants[h_key]
        x_values = np.arange(len(h_data))
        plt.close('all')  # Close all existing figures
        fig, ax = plt.subplots()
        markerline, stemlines, baseline = ax.stem(x_values, h_data, '-.')
        plt.setp(baseline, 'color', 'k', 'linewidth', 2)
        plt.title('Stem plot of selected filter')
        plt.xlabel('Index')
        plt.ylabel('Amplitude')
        plt.show()

# Setup the widgets
dropdown = widgets.Dropdown(options=list(h_variants.keys()), value='h32', description='Filter:')

# Function that updates the plot based on dropdown
def update_plot(change):
    plot_stem(change['new'])

# Observe dropdown for changes
dropdown.observe(update_plot, names='value')

# Group the dropdown and output widget in a vertical box layout
vbox = widgets.VBox([dropdown, output2])

# Display the VBox
display(vbox)

# Initial call to display the plot
plot_stem(dropdown.value)  # Use the dropdown's value
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
import ipywidgets as widgets
from IPython.display import display

# Compute frequency responses
freq32, resp32 = signal.freqz(h32)
freq64, resp64 = signal.freqz(h64)
freq128, resp128 = signal.freqz(h128)
freq160, resp160 = signal.freqz(h160)
freq256, resp256 = signal.freqz(h256)

# Function to apply window and compute frequency responses
def compute_filtered_freq_responses(window_type='Rectangular'):
    freqs, resps = {}, {}
    for filt_size in [32, 64, 128, 160, 256]:
        filt_name = f'h{filt_size}'
        filt = eval(filt_name)
        
        if window_type == 'Hamming':
            filt = filt * np.hamming(filt_size)
        elif window_type == 'Kaiser':
            beta = 14  # Example beta value for Kaiser window, adjust as needed
            filt = filt * np.kaiser(filt_size, beta)
            
        freqs[filt_name], resps[filt_name] = signal.freqz(filt)
    return freqs, resps

# Initial computation with Rectangular (no window)
freqs, resps = compute_filtered_freq_responses()

output3 = widgets.Output()

def update_plot(change=None):
    window_type = window_type_dropdown.value
    freqs, resps = compute_filtered_freq_responses(window_type)
    
    with output3:
        clear_output(wait=True)
        plt.figure(figsize=(10, 5))
        for filt in ['h32', 'h64', 'h128', 'h160', 'h256']:
            if checkboxes[filt].value:
                w, h = freqs[filt], resps[filt]
                plt.plot(0.5 * Fs * w / np.pi, 20 * np.log10(np.abs(h)), label=filt)
        plt.title(f'Frequency Response with {window_type} Window')
        plt.xlabel('Frequency (Hz)')
        plt.ylabel('Magnitude (dB)')
        plt.xscale('linear')
        plt.grid(True)
        plt.legend()
        plt.show()

# Create checkboxes for each filter
checkboxes = {f'h{size}': widgets.Checkbox(value=False, description=f'h{size}') for size in [32, 64, 128, 160, 256]}
for cb in checkboxes.values():
    cb.observe(update_plot, names='value')

# Dropdown for selecting the window type
window_type_dropdown = widgets.Dropdown(
    options=['Rectangular', 'Hamming'],
    value='Rectangular',
    description='Window Type:',
)
window_type_dropdown.observe(update_plot, names='value')

# VBox to hold the checkboxes and dropdown
vbox_checkboxes = widgets.VBox(list(checkboxes.values()) + [window_type_dropdown])

# Display the VBox along with the output area
vbox = widgets.VBox([vbox_checkboxes, output3])
display(vbox)
  1. Σχεδιάστε φίλτρο Parks-McClellan μήκους 160+1, με τις ίδιες οριακές συχνότητες, όπως στο παράδειγμα (0.1, 0.15). Σχεδιάστε την απόκριση συχνότητας και αυτού του φίλτρου και σχολιάστε τις διαφορές από το ίδιο φίλτρο μήκους 64+1 (του παραδείγματος 1.2). Γράψτε την απάντησή σας στο αρχείο κειμένου lab2_nnnnn.txt.

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
import ipywidgets as widgets
from IPython.display import display, clear_output

filters = {
    'Equiripple Filter 32+1': signal.remez(33, [0, 0.1*Fs, 0.15*Fs, 0.5*Fs], [1, 0], Hz=Fs),
    'Equiripple Filter 64+1': signal.remez(65, [0, 0.1*Fs, 0.15*Fs, 0.5*Fs], [1, 0], Hz=Fs),
    'Equiripple Filter 128+1': signal.remez(129, [0, 0.1*Fs, 0.15*Fs, 0.5*Fs], [1, 0], Hz=Fs),
    'Equiripple Filter 160+1': signal.remez(161, [0, 0.1*Fs, 0.15*Fs, 0.5*Fs], [1, 0], Hz=Fs),
    'Equiripple Filter 256+1': signal.remez(257, [0, 0.1*Fs, 0.15*Fs, 0.5*Fs], [1, 0], Hz=Fs)
}

output_plot = widgets.Output()

# Function to update the plot based on selected filters
def update_plot(change):
    selected_filters = [cb.description for cb in checkboxes if cb.value]
    with output_plot:
        clear_output(wait=True)
        plt.figure(figsize=(10, 5))
        for filter_name in selected_filters:
            filter_coeffs = filters[filter_name]
            w, h = signal.freqz(filter_coeffs, worN=8000)
            plt.semilogy(w * Fs / (2 * np.pi), np.abs(h), label=filter_name)
        plt.title('Frequency Response')
        plt.xlabel('Frequency (Hz)')
        plt.ylabel('Gain')
        plt.legend()
        plt.grid(True)
        plt.show()

# Create checkboxes for each filter
checkboxes = [widgets.Checkbox(value=False, description=name) for name in filters.keys()]

# Register the update function with each checkbox
for cb in checkboxes:
    cb.observe(update_plot, names='value')

# VBox to hold the checkboxes and output
vbox = widgets.VBox([widgets.Label('Select Equiripple Filter Length:')] + checkboxes + [output_plot])

# Display the VBox
display(vbox)
  1. Αλλάξτε τις οριακές συχνότητες του φίλτρου του ερωτήματος 3 σε (0.11, 0.12) και συγκρίνετε τις αποκρίσεις συχνότητας των δύο φίλτρων. Αποθηκεύστε τα σχεδιαγράμματα και γράψτε τα σχόλιά σας στο αρχείο κειμένου lab2_nnnnn.txt.

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
import ipywidgets as widgets
from IPython.display import display, clear_output

filters = {
    'Equiripple Filter 32+1': signal.remez(33, [0, 0.1*Fs, 0.15*Fs, 0.5*Fs], [1, 0], Hz=Fs),
    'Equiripple Filter 64+1': signal.remez(65, [0, 0.1*Fs, 0.15*Fs, 0.5*Fs], [1, 0], Hz=Fs),
    'Equiripple Filter 128+1': signal.remez(129, [0, 0.1*Fs, 0.15*Fs, 0.5*Fs], [1, 0], Hz=Fs),
    'Equiripple Filter 160+1': signal.remez(161, [0, 0.1*Fs, 0.15*Fs, 0.5*Fs], [1, 0], Hz=Fs),
    'Equiripple Filter 256+1': signal.remez(257, [0, 0.1*Fs, 0.15*Fs, 0.5*Fs], [1, 0], Hz=Fs)
}

output_plot = widgets.Output()

# Function to update the plot based on selected filters
def update_plot(change):
    selected_filters = [cb.description for cb in checkboxes if cb.value]
    with output_plot:
        clear_output(wait=True)
        plt.figure(figsize=(10, 5))
        for filter_name in selected_filters:
            filter_coeffs = filters[filter_name]
            w, h = signal.freqz(filter_coeffs, worN=8000)
            plt.semilogy(w * Fs / (2 * np.pi), np.abs(h), label=filter_name)
        plt.title('Frequency Response')
        plt.xlabel('Frequency (Hz)')
        plt.ylabel('Gain')
        plt.legend()
        plt.grid(True)
        plt.show()

# Create checkboxes for each filter
checkboxes = [widgets.Checkbox(value=False, description=name) for name in filters.keys()]

# Register the update function with each checkbox
for cb in checkboxes:
    cb.observe(update_plot, names='value')

# VBox to hold the checkboxes and output
vbox = widgets.VBox([widgets.Label('Select Equiripple Filter Length:')] + checkboxes + [output_plot])

# Display the VBox
display(vbox)

Εφαρμογή του φίλτρου#

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.io.wavfile import write
from IPython.display import display, Audio, clear_output
import ipywidgets as widgets
import io

lp_filters = {
    'lpass32': signal.remez(33, [0, 0.1*Fs, 0.15*Fs, 0.5*Fs], [1, 0], Hz=Fs),
    'lpass64': signal.remez(65, [0, 0.1*Fs, 0.15*Fs, 0.5*Fs], [1, 0], Hz=Fs),
    'lpass128': signal.remez(129, [0, 0.1*Fs, 0.15*Fs, 0.5*Fs], [1, 0], Hz=Fs),
    'lpass160': signal.remez(161, [0, 0.1*Fs, 0.15*Fs, 0.5*Fs], [1, 0], Hz=Fs),
    'lpass256': signal.remez(257, [0, 0.1*Fs, 0.15*Fs, 0.5*Fs], [1, 0], Hz=Fs)
}

# Prepare audio data
audio_data = {}

for name, lp_filter_coeffs in lp_filters.items():
    # Apply the filter to the signal
    s_lp_filtered = signal.convolve(s, lp_filter_coeffs, mode='same') / np.sum(lp_filter_coeffs)
    
    # Normalize the filtered signal to prevent clipping
    s_lp_filtered_normalized = np.int16((s_lp_filtered / np.max(np.abs(s_lp_filtered))) * 32767)
    
    # Write to an in-memory file
    audio_buf = io.BytesIO()
    write(audio_buf, Fs, s_lp_filtered_normalized)
    audio_buf.seek(0)  # Go back to the beginning of the BytesIO object
    
    # Store the audio buffer
    audio_data[name] = audio_buf

# Function to display the frequency response of the filtered signal
def plot_freq_response(filter_name):
    with plot_output:
        clear_output(wait=True)
        # Apply the filter to the original signal
        s_lp_filtered = signal.convolve(s, lp_filters[filter_name], mode='same') / np.sum(lp_filters[filter_name])
        
        # Calculate FFT
        freqs = np.fft.rfftfreq(len(s_lp_filtered), 1/Fs)
        fft_mag = np.abs(np.fft.rfft(s_lp_filtered))
        
        # Convert magnitude to dB
        fft_mag_db = 20 * np.log10(fft_mag)
        
        plt.figure(figsize=(10, 5))
        plt.plot(freqs, fft_mag_db)
        plt.title(f'Frequency Response of Filtered Signal with {filter_name}')
        plt.xlabel('Frequency (Hz)')
        plt.ylabel('Magnitude (dB)')
        plt.grid(True)
        plt.xlim(0, Fs/2)
        plt.show()

# Audio playback function
def show_audio_player(filter_name):
    with audio_output:
        clear_output(wait=True)
        display(Audio(data=audio_data[filter_name].getvalue(), rate=Fs))

# Widget setup
filter_dropdown = widgets.Dropdown(
    options=[(name, name) for name in lp_filters.keys()],
    value='lpass32',
    description='Filter:'
)

plot_output = widgets.Output()
audio_output = widgets.Output()

def on_filter_change(change):
    filter_name = change['new']
    plot_freq_response(filter_name)
    show_audio_player(filter_name)

filter_dropdown.observe(on_filter_change, names='value')

# Display the UI
vbox = widgets.VBox([filter_dropdown, plot_output, audio_output])
display(vbox)

# Initialize
on_filter_change({'new': filter_dropdown.value})
  1. Αντικαταστήστε το σήμα s με άθροισμα τεσσάρων ημιτονικών συναρτήσεων μοναδιαίου πλάτους, συχνότητας 700, 900, 1400 και 2500 Hz και διάρκειας 1.0 sec. Διατηρήστε την ίδια Fs = 8192 Hz. Σχεδιάστε τη φασματική πυκνότητα του σήματος (με pwelch) και δείτε το αποτέλεσμα του φιλτραρίσματος με τα δύο φίλτρα Parks-McClellan των ερωτημάτων 3 και 4. Αποθηκεύστε τα φασματικά διαγράμματα στο αρχείο κειμένου lab2_nnnnn.txt.

Μέρος 3: Εφαρμογή A#

Εκτελέστε την άσκηση 1.2 της παραγράφου 1.6 του τεύχους του μαθήματος, μεταβάλλοντας κατάλληλα τον κώδικα του παραδείγματος 1.2 (που έχετε αποθηκεύσει στο αρχείο lab2_1_nnnnn.m). Σκοπός της άσκησης είναι ο σχεδιασμός ενός ζωνοπερατού φίλτρου ζώνης διέλευσης (0.7 ΚHz, 1.5 ΚHz) με τις δύο μεθόδους που περιγράφηκαν πιο πάνω και η εφαρμογή τους στο σήμα s του παραδείγματος 1.2 (δείτε και το Παράδειγμα 2 στο τέλος της παρουσίασης «Η ΨΕΣ στις τηλεπικοινωνίες»). Επαληθεύστε το σωστό σχεδιασμό του φίλτρου σας, ελέγχοντας τόσο την απόκριση συχνότητάς του όσο και το αποτέλεσμα του φιλτραρίσματος στο σήμα s

Ζωνοπερατά φίλτρα#

Με αναλυτικό υπολογισμό της κρουστικής απόκρισης και παράθυρο#

Hide code cell source
# Με αναλυτικό υπολογισμό της κρουστικής απόκρισης και παράθυρο kaiser
f1=700; f2=1500;  
Ts=1/Fs
f2m1=(f2-f1); f2p1=(f2+f1)/2; N=256
t=np.arange(-(N-1),N-1,2)*Ts/2
hbp=2/Fs*np.divide(np.multiply(np.cos(2*np.pi*f2p1*t),np.sin(np.pi*f2m1*t))/np.pi,t)
hbpw=np.multiply(hbp,signal.kaiser(len(hbp),5))

s1_bp=signal.convolve(s1,hbp,'same')
Hide code cell source
f0, Pxx1_den = signal.welch(s1_bp, Fs, noverlap=128, nperseg=256)

# Create Plotly figure
fig = go.Figure()

# Add trace for the power spectral density
fig.add_trace(go.Scatter(x=f0, y=Pxx1_den, mode='lines', line=dict(color='blue')))

# Update layout for a semilog plot
fig.update_layout(
    title='Φάσμα ζωνοπερατού σήματος sonar',title_x=0.5,title_font=dict(size=20, color='black', family="Arial, sans-serif"),
    xaxis_title='Συχνότητα (Hz)',
    yaxis_title='Πυκνότητα φάσματος ισχύος',
    yaxis_type='log',  # Set y-axis to logarithmic scale
    template='plotly_white',
    xaxis=dict(showgrid=True),
    yaxis=dict(showgrid=True, range=[np.log10(1e-15), np.log10(1e-7)])  # Adjust y-axis range
)

# Show figure
fig.show(renderer='notebook')


def play_sound(b):
    sd.play(20 * s1_bp, Fs)

# Create a button widget
play_button = widgets.Button(description="Play Sound")

# Link the button to the function that plays the sound
play_button.on_click(play_sound)

# Display the button
display(play_button)

Ζωνοπερατό ισουψών κυματώσεων#

Hide code cell source
bpass = signal.remez(128, [0, f1*0.95, f1*1.1, f2*0.95, f2*1.05, Fs/2], [0, 1, 0], fs=Fs)
freq,resp_pm = signal.freqz(bpass)

# Compute frequency response for the equiripple bandpass filter
freq, resp_pm = signal.freqz(bpass)

# Create Plotly figure
fig = go.Figure()

# Add trace for equiripple bandpass filter response
fig.add_trace(go.Scatter(x=0.5 * Fs * freq / np.pi, y=np.abs(resp_pm), mode='lines', name='Equiripple Bandpass', line=dict(color='green')))

# Update layout for a semilog plot
fig.update_layout(
    title='Απόκριση συχνότητας ζωνοπερατού φίλτρου equirriple',title_x=0.5,title_font=dict(size=20, color='black', family="Arial, sans-serif"),
    xaxis_title='Συχνότητα (Hz)',
    yaxis_title='Κέρδος',
    yaxis_type='log',  # Set y-axis to logarithmic scale
    template='plotly_white',
    xaxis=dict(showgrid=True),
    yaxis=dict(showgrid=True)
)

# Show figure
fig.show(renderer='notebook')

# Convert figure to HTML
fig_html = pio.to_html(fig, full_html=False, include_plotlyjs='cdn')
display(HTML(fig_html))
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[15], line 28
     25 fig.show(renderer='notebook')
     27 # Convert figure to HTML
---> 28 fig_html = pio.to_html(fig, full_html=False, include_plotlyjs='cdn')
     29 display(HTML(fig_html))

NameError: name 'pio' is not defined

Εδώ θα πειραματιστούμε με δύο σήματα:

(i) το sonar του παραδείγματος, το οποίο εδώ διαβάζεται από ένα .txt αρχείο (έχει προέλθει με εξαγωγή του s από το MATLAB) και

Στο παράρτημα:

(ii) ένα σήμα μουσικής, το violin.wav (σήμα από μουσική βιολιού), το οποίο περιέχει υψηλότερες συχνότητες και έχει προέλθει με δειγματοληψία στα Fs_viol=44100 Hz.